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Abstract. Comparison is made between a number of independent computer programs for radiative transfer in 
molecular rotational lines. The test models are spherically symmetric circumstellar envelopes with a given density 

' and temperature profile. The first two test models have a simple power law density distribution, constant temper- 

ature and a Active 2-level molecule, while the other two test models consist of an inside-out collapsing envelope 
observed in rotational transitions of HCO + . For the 2-level molecule test problems all codes agree well to within 
0.2%, comparable to the accuracy of the individual codes, for low optical depth and up to 2% for high optical 
depths (r=4800). The problem of the collapsing cloud in HCO + has a larger spread in results, ranging up to 
12% for the J— A population. The spread is largest at the radius where the transition from collisional to radiative 
excitation occurs. The resulting line profiles for the HCO + J=4-3 transition agree to within 10%, i.e., within the 
• »"H . calibration accuracy of most current telescopes. The comparison project and the results described in this paper 

' provide a benchmark for future code development, and give an indication of the typical accuracy of present day 

calculations of molecular line transfer. 
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1. Introduction 

Molecular lines are excellent probes of the physical and 
chemical conditions in interstellar clouds, protostellar 
envelopes, circumstellar shells around late-type stars, 
photon-dominated regions etc. Furthermore, molecular 
line transitions play a key role in probing the properties 
of galaxies and their evolution. The interpretation of such 
lines requires the use of line radiative transfer programs 
which can calculate accurately the non-LTE (local thermo- 
dynamic equilibrium) level populations and the resulting 
output spectra. See Black (2000) for a recent review. 

It is known from stellar atmosphere research that sub- 
tle errors in radiative transfer algorithms can lead to sig- 
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nificantly incorrect results (Mihalas 1978). A particularly 
well-known problem of this kind is an insufficiently strin- 
gent convergence criterion at high optical depths. In the 
absence of a-posteriori checks on numerical results, the 
best way to validate the methods is by the use of vari- 
ous techniques, if possible with many independent codes 
of different types. Once the reliability and the behavior of 
the codes has been established and understood, they can 
be safely used within the limits of parameter space within 
which tests have been carried out. 

In this paper, the setup and results of a large-scale 
campaign to compare line radiative transfer programs are 
described. The test problems are spherically symmetric, 
and can be modeled by a 1-dimensional radiative transfer 
code. Codes capable of handling more than one dimension 
could be compared in a similar way in the future. 

The main aim of this paper is to examine and com- 
pare the various radiative transfer methods that are cur- 
rently used in the astrophysical community for modeling 
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the shape and strength of molecular lines emerging from 
interstellar and circumstellar matter. A series of prob- 
lems is chosen that represent the difficulties that may 
be encountered in comparison with data from present- 
day and future ground-based telescopes, such as the Sub- 
Millimeter Array (SMA), the Atacama Large Millimeter 
Array (ALMA) and future airborne and space-borne mis- 
sions such as SOFIA and the Herschel Space Observatory. 
The test problems and their solutions presented here are 
available to the community via a web-site, allowing users 
to run the same problems, check the accuracy of their 
codes and thus speed up the further development of their 
own radiative transfer codes. It is hoped that this will 
stimulate a more widespread use of these codes for the 
interpretation of molecular line observations. 

We present two test problems, each at two differ- 
ent optical depths (i.e. in total four problems). The first 
problem is a fictive 2-level molecule in a spherical enve- 
lope with a powerlaw density with no systematic veloc- 
ity and a constant temperature. This problem is meant 
to tune all codes before examining the more complex, 
realistic problem. The main test problem is based on 
the 1-dimensional (1-D) inside-out collapse model, where 
the level populations are computed for HCO + at vari- 
ous abundances. The HCO + ion is chosen as a represen- 
tative example of the molecule which samples gas with 
a large range in densities and is readily observable in a 
variety of astrophysical regions. Both a high and a low 
optical depth model, representative of the main isotope 
HCO + and the less abundant isotope H 13 CO + , are used 
to check the convergence properties of each code. All test 
problems and their results can be found at the webpage 
http : //www . strw . leidemmiv . nl/ ^radtrans . 



2. Molecular line radiative transfer 

The radiative transfer problem is represented by an equa- 
tion describing the emission, absorption and movement of 
photons along a straight line in a medium: 

= 2 v -a v I v . (1) 



with 



the 



ds 
notation 



adopted 

Rybicki fc Lightman (1979)| . Equation (§) is 



from 
the dif- 
ferential description of the intensity I v along a photon 
path (ds) at frequency v, where j„ [erg s -1 cm -3 ] and a v 
[cm -1 ] denote the emission and absorption coefficients. 
Another way of writing Eq. ([!]) is: 



dl^ q _ r 

GIT,, 



(2) 



where S v — ju/oi v is referred to as the source function 
(the emissivity of the medium per unit optical depth), 
and the optical depth r v is defined in differential form as 
g?t„ = a v ds. Equation (||) can be written in integral form, 
which is the form that is most often used in radiative 
transfer codes: 



S„{T')e T '- T dr' . 



(3) 



where r is the optical depth between the point where I v is 
evaluated and spatial infinity along the line (i.e., s — — oo). 
This integral is evaluated along all possible straight lines 
through the medium. In practice, these will be a discrete 
sample of lines covering space and direction as well as 
possible. 

For the problem of molecular line transfer the emission 
and absorption coefficients are determined by the transi- 
tion rates between the various rotational and/or vibra- 
tional levels of the molecule, and the population of these 
levels. For a transition from level i to level j (where the 
energy of level i is greater than that of level j) the emission 
and absorption coefficients are given by: 



jij{v) = niAij4>ij(v) , 



(4) 
(5) 



where rii and rij [cm -3 ] are the population densities of the 
upper (i) and lower (j) level, and Aij, Bij and Bji are the 
Einstein coefficients. The function <fiij(v) represents the 
line profile for this transition, which is properly described 
as a Voigt profile, a combination of a (micro-turbulent) 
Gaussian and intrinsic Lorentzian line broadening. 

The source function for a particular transition Sij is 
independent of velocity if one assumes complete frequency 
redistribution, i.e., the frequency deviation from line cen- 
ter of absorbed and emitted photons are uncorrelated. The 
source function then becomes: 



lJ aij(v) UjBji-riiBi 



(6) 



The relative level populations n-i are determined from 
the statistical equilibrium equation: 



^[rijAji + (rijBji - mB^Jji] 



3>l 
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(7) 



where C%j = n co \Kij with K%j the collisional rate coeffi- 
cients in cm 3 s _1 and n co \ the density of collision partners, 
taken here to be H2 in J=0. Jji is the integrated mean 
intensity over the line profile: 



Jij = — I I U (Q) <t>{y) dfl dv . 



(8) 



The symbol represents the spatial direction in which 
the intensity I V (Q) is measured. 

A useful concept for the analysis of radiative trans- 
fer calculations is the excitation temperature T ex of the 
transition between level i and level j, given by 



— Kva. 



in [ am 

9i n 3 



(9) 



where k is Boltzmann's constant and Qi the statistical 
weight of level i. The energy difference between the two 
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Eq 3: Radiative transfer equation 



Sv 



Eq 8: Mean intensity calculation 



Eq 5: oty calculation 

Eq 6: Source function calculation 



Eq 7: Statistical equilibrium eq. 



Fig. 1. Flow diagram of the molecular line radiative trans- 
fer problem. Since the problem is coupled (follow arrows), 
several iterations are needed to calculate the true level 
populations rii. All adopted symbols are explained in §|| 

levels is given by hvij . In local thermodynamic equilibrium 
(LTE) , T cx equals the local gas temperature, while if T ex 
higher or lower, the excitation is super- or sub-thermal. In 
addition, the intensity is proportional to T ox in the opti- 
cally thin limit. 

Equations form a coupled set of equations. The 

way in which the quantities depend on each other is de- 
picted in Fig. [jj. The intensities are found by integrating 
the source function (Eq. 3). The source function and ex- 
tinction coefficients depend on the level populations (Eqs. 
4-5). These in turn depend on the intensities (Eqs. 7- 
8). To solve this set of equations one must determine the 
radiation field and the level populations simultaneously. 
Since the radiation field couples the level populations at 
different spatial positions to each other through the trans- 
fer integral (Eq. ||), the only way of solving this system 
directly seems to be through complete linearization and 
solving a huge matrix equation involving all level popu- 
lations at all spatial positions. In practice, the evaluation 
of the matrix elements and the inversion of the problem 
consumes far too much CPU-time as well as memory and 
is therefore beyond current computing capabilities. An al- 
ternative and much simpler way is to iteratively evalu- 
ate all equations, following the arrows, as illustrated in 
Fig. [lj. This method is called "Lambda Iteration" and in- 
cludes formal integral codes as well as Monte Carlo meth- 
ods. Though simple, its principal disadvantage is that it 
converges slowly at high optical depths. Many radiative 
transfer codes therefore use a hybrid scheme: the direct 
inversion of a simplified subset of the equations, and itera- 
tion to solve the remaining problem. In this paper we have 
used codes of the Lambda Iteration type (including Monte 
Carlo methods), of the "Complete Linearization" type, 
and hybrid schemes such as "Approximate/ Accelerated 
Lambda Iteration" and "Accelerated Monte Carlo" . 

3. Methods and codes 

In this paper four methods and eight different codes are 
compared. Any practical method discretizes the problem 
to determine the level populations at each position. A 



short description is given in this section of each of these 
codes. The first method (LI) is not used, but given as an 
introduction for the other methods. In Table fll the dif- 
ferent codes and principal authors are indicated. Detailed 
descriptions of each of the codes can be found in the ref- 
erences given. 

3.1. Lambda Iteration (LI) 

The Lambda Iteration method is the most basic method 
of all. It involves the iterative evaluation of level popula- 
tions and intensities until the system has converged. The 
name "Lambda Iteration" originates from the fact that the 
process of iteration can be mathematically written into a 
formalism involving a "Lambda Operator" . This Lambda 
Operator represents the entire procedure of computing the 
line-weighted mean intensities Jy from the source func- 
tion. It involves the formal integral Eq. (||) along all pos- 
sible lines through the medium, and includes the angle- 
frequency integrations of Eq. (||) to obtain the J. The 
Lambda Operator is defined as: 



J ij A-ij [Sij ] , 



(10) 



and is therefore a global operator. The A„ operator is a 
matrix operator connecting all points and all levels to 
each other. To solve the radiative transfer numerically, 
the problem has to be discretized both in space and fre- 
quency and initial level populations have to be assumed. 
With these level populations, the Lambda Operator is 
constructed by solving the radiative transfer for I v and 
J, which is inserted in Eq. (0) to calculate the new level 
populations. The procedure is repeated until the relative 
change in the level populations or mean intensities be- 
tween two successive steps falls below a desired conver- 
gence criterium. The time to reach convergence of the so- 
lution is proportional to r 2 for r ^ 1 and is therefore ex- 
tremely slow for highly optically thick regions. The very 
small local changes in the level populations in successive 
iterations could then be easily mistaken for convergence. 

3.2. Monte Carlo (MC) 

The Monte Carlo (MC) method is based on the simula- 
tion of basic physical processes with the aid of random 
numbers. This makes the formalisms of these codes rela- 
tively simple and intuitive as one has to deal only with 
the basic formulae. But one must take proper care of 
the statistics and be sure that all regions in space and 
frequency are well sampled by the randomly distributed 
photon packages. The Monte Carlo method for molecu- 
lar line transfer has been described by Berncs (1979)| for 
a spherically symmetric cloud with a uniform density. A 
major advantage of MC codes is the possibility of non- 
regular grids, in particular for multi-dimensional prob- 
lems. In fact, MC methods are straightforward to extend 
from 1-D to 2-D (Hogerheijde & van der Tak 2000) and 
even 3-D (Park & Hong 1998, Juvela 1997). One of the 
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major disadvantages of the method is the long CPU time 
which is needed to lower the random error intrinsic in 
the method. This error decreases inversely proportional 
to the square root of the number of simulated photons. 
In addition, the method suffers from similar convergence 
problems as Lambda Iteration, as it is formally a variant 
of the Lambda Iteration method. 

Thus far, the Monte Carlo method has been imple- 
mented in two ways. In one set of codes, the radiation 
field is randomly sampled by discrete photons that are 
emitted and absorbed in the material. In the other set, 
the radiation field is sampled by random directions along 
which the radiative transfer is evaluated. 

3.3. Approximated/ Accelerated Lambda Iteration 
(All) 

This scheme is similar to the Lambda iteration scheme, 
with the difference that the equations are pre-conditioned 
to speed up the convergence. The simplest (but still highly 
efficient) way to implement the ALI scheme is by splitting 
the Lambda Operator into the local self-coupling contri- 
bution and the remainder: 



A lJ = (A 4 ,-A* 7 ) + A* 



(11) 



with Ay the lambda operator between the i and j levels 
and A*j the diagonal (local), or tridiagonal (local -l-nearest- 
neighbor), part of the full Lambda operator. In general 
more reliable convergence properties are obtained with 
tridiagonal or higher bandwidths (Hauschildt et al. 1994). 
Since an approximate operator A* 3 of this kind is easily 
invertible, and since its matrix elements are relatively easy 
to calculate, one can use this operator for the matrix in- 
version, and use the iteration for solving the remainder 
(non-local part) of the transfer problem. Since the self- 
coupling of a cell at high optical depth is the bottle-neck 
that slows down convergence in Lambda Iteration codes, 
this removal of the local self-coupling by direct inversion 
of A*j can speed up the convergence drastically. A full 
description of the method was given by, e.g., Rybicki & 
Hummer (1991, 1992). 

In addition to this operator splitting technique, one 
can apply certain iteration-improvement schemes such as 



the Ng-acceleration method (Ng 1974). These are methods 
that can improve the convergence of any linearly converg- 
ing iteration scheme. In the Ng-method one uses the pre- 
vious results to estimate the convergence behavior of the 
problem. After every four iteration steps an extrapolation 
can then be performed towards the expected convergence. 
The number of iterations is not a strict requirement but 
four is typically found to give a reliable and significant 
acceleration. This scheme is very effective as can be seen 
in Figs. 1 to 3 in Rybicki & Hummer (1991), where the 
number of iterations is plotted versus the convergence. 



3.4. Accelerated Monte Carlo (AMC) 

The difference in Accelerated Monte Carlo and Monte 
Carlo methods can be similarly described as the ALI com- 
pared to LI method. Most important is the splitting of the 
mean intensity in an external field and a local contribu- 
tion. Formally, 



J = (A-A*)[^(J)] + A*[S ui (J)] 



—external —local 

= J +J 



N ^ 



(12) 
(13) 
(14) 



with S^J) the results from the previous iteration and N 
the number of photons. This is an important issue as in 
the standard Monte Carlo approach, most time is taken 
by photons trapped in an optically thick cell, where due 
to absorption and randomly directed emission the photons 
follow a random walk path through the cell, with one step 
for each iteration. 

The splitting can be done in different ways. Juvcla 
(1997) implements a diagonal lambda operator in Eqs. 
(12)-(14) by counting explicitly those photons which were 
emitted and absorbed within the same cell. A reference 
field described by Bernes et al. (1979) is used to decrease 
the random fluctuations caused by the Monte Carlo sam- 
pling. Hogerheijde & van der Tak (2000) perform sub- 
iterations to calculate the local contribution in each cell, 
thus ensuring that in each cell the local field and popula- 
tions are always consistent with the external field due to 
all other cells. A third possibility of acceleration, adopted 
by Schoier (2000), is the use of core-saturation (Rybicki 
1972, Hartstein & Liseau 1998), where optically thick line 
center photons are replaced by the local source function. 

3.5. Local Linearization (MULTI type codes) 

This method was developed by 

|Scharmer fc Carlsson (1985) and Harper (1994) to 
produce the MULTI and SMULTI line radiation trans- 
port codes respectively. This approach perturbs Eqs. (2) 
to (8) linearly, neglecting second order or higher terms. 
The linearization greatly reduces the effect of high optical 
depth terms and allows rapid convergence. The MULTI 
and SMULTI codes use the Olson diagonal approximate 
operator (Olson, Aucr, & Buchlcr 1986) scheme to save 
on storage. This operator uses the diagonal of the A 
matrix when computing the solution to the linearized 
equations. Although this adds in approximations and 
delays convergence, it greatly reduces the required storage 
capacity and so for many problems it allows the problem 
to be tackled in the first place. Olson's scheme simply 
assumes that the changes in the intensity are related to 
the changes in the source and opacity terms (the level 
populations). The main difference to the other schemes 
is that the solution to the linearized equations returns 
changes, Sn and 5 J to the level populations and the mean 
intensities respectively. There is always the possibility 
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that a converged solution contains unphysical negative 
populations, and this provides a useful indicator for 
poor sampling and errors. In the version used here a 
ALI-type scheme to complement the MULTI is used, 
providing numerical stability but also slower convergence 
in optically thick media. 

3.6. Radiative transfer codes 



Each code described in this paper uses one or more of 
the above techniques to calculate a converged set of level 
populations. In this section we list the codes which use 
each of the convergence methods listed above. Then for 
each code we describe how convergence is accelerated, how 
each code samples the volume under consideration, and 
state the convergence criteria used by each code. 

1. Monte Carlo ( MC): 

— F. Schoier ( ^choier 2000| ); The rate of convergence 
depends on the number of model photons and 
the iterative procedure. In the problems presented 
here, the counters of stimulated emission are reset 
after each set of 5 iterations. The number of it- 
erations needed for convergence in the "classical" 
Monte Carlo scheme is of the same order as the 
maximum optical depth in the model. The core 
saturation method is included to speed up the con- 
vergence in the high optical depth case. To reach 
an accuracy of ^10 -2 in the derived level popula- 
tions, ~ 10 5 — 10 6 model photons per iteration are 
generally needed in the optically thick case. 

2. Accelerated Lam bda Iteration (ALI): 

— V. Ossenkopf flOsscnkopf, Trojan, fc Stutzki 2001 ); 
All discretizations are done by the code such that 
sufficiently small differences between the grid 
points are guaranteed. This concerns the radial 
shells, the grid of rays, the frequency grid, and 
the number of levels used. The convergence is 
measured in terms of radiative energy densities, 
not level populations, which are extrapolated 
from the Auer acceleration scheme. Locally, 
the code may take into account turbulent sub- 
clumping using the statistical description of 
Martin, Sanders fc Hills (1984)| , but this is not 
used here. 

— S. Doty; the input data are interpolated onto the 
spatial grid of the problem with a variable method 
(log/linear/combined). This, coupled with extreme 
care in the ray-tracing integration, helps to ensure 
high accuracy. A local approximate lambda opera- 
tor and Ng acceleration are used to solve both the 
line and continuum transfer. Convergence is deter- 
mined by changes in both level populations and the 
local net heating rate, with small (< 10~ 4 ) changes 
required for the entire spatial grid for both accel- 
erated and non-accelerated iterations. 



— CP. Dullemond ( pullcmond fc Turolla 2000| ); 
The formal integration proceeds via a short- 
characteristics method using the three-point 
quadrature formula of Olson & Kunasz (1987). 
The discrete local angles fii are chosen using 
the tangent-ray prescription. This means that 
successive short characteristics match up and no 
interpolations are required. At every grid point 
there are 41 angular points in fi with three extra 
around [i ~ to prevent undersampling. Using 
the ALI method with a local operator and Ng 
acceleration, the system is iterated towards a 
solution with a convergence criterion 10 -6 in level 
population. 

— H. Wiesemeyer ( |Wicsemeyer 1997) The rate of 
convergence is improved by applying the method of 
minimization of residuals (Auer 1987). Other accel- 
eration methods, such as vector and matrix extrap- 
olation, are currently being implemented. The so- 
lution to the equation of radiative transfer is evalu- 
ated by various numerical methods (either quadra- 
ture rules, or multi- value methods), according to 
the accuracy required by the problem to be solved. 
The code uses long characteristics, to preserve an 
isotropic distribution of rays at any point of the 
model volume, following the multivariate quadra- 
ture rule of Steinacker et al. (1996). Its performance 
is thus rather optimized to solve smooth problems 
in more than one dimension. All discretizations are 
performed by the code. 

Accelerated Monte Carlo (AMC): 

— M. Juvela (Juvela 1997); This code uses the refer- 
ence field method to reduce random fluctuations. 
The velocity is discretized into 50 channels and the 
number of rays per iteration was taken to be 350. 
The ray generation is weighted such that the same 
number of rays was shot through each annulus. The 
random number generators are reset after each it- 
eration, making it possible to use Ng acceleration 
on every third iteration starting with the fifth it- 
eration. Convergence was tested only for the eight 
lowest levels down to lxlO -3 in level population. 

— M. Hogcrhcijde fc F. v an der Tak 
; van der Tak 200C); In this code, 

time is considered, with TV rays 
from random directions. The 



(rlogerheijde 
one cell at a 
entering the cell 
radiative transfer is followed along each of these 
rays, starting with the CMB field at the edge of 
the cloud. In this way it is possible to calculate 
separately the local contribution to the radiation 
field in the target cell, allowing significant re- 
duction of the computing time for optically thick 
cells. For a first order estimate of the radiation 
field, the same set of random numbers is used 
to describe the ray directions for each iteration, 
thus resembling a fixed-ray (ALI) code. This 
process yields a solution free of random noise but 
with possible inadequate sampling of directions 
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Table 1. Codes used 



Label 


Author 


method 


dimensions 


Keicrenec 


A 


Juvela 


AMC 


ID, 2D & 3D 


|. Juvela (1997) 


B 


Hogerheijde & 


AMC 


ID & 


i 2D 


Hogerheijde & van der Tak (2000) 




van der Tak 










C 


Schoier 


MC 


ID 




Schoier (2000) 


D 


Doty 


ALI 


ID 






E 


Ossenkopf 


ALI 


ID & 


i 2D 


Ossenkopf, Trojan & Stutzki (2001) 


F 


Dullemond 


ALI 


ID & 2D 


Dullemond & Turolla (2000) 


G 


Yates 


(S)MULTI 


ID 




Rawlings & Yates (2001) 


H 


Wiesemeyer 


ALI 


ID & 2D 


Wiesemeyer (1997) 



1 Label used in Figures 4 and 6 for each of the codes. 



and velocities. In the second stage a different set 
of random numbers is used for each iteration, 
providing true random sampling of the radiation 
field. The number of rays is increased for each 
cell until convergence is reached, ensuring proper 
angular sampling of the radiation field everywhere. 
The equations of statistical equilibrium are solved 
in each iteration to a fractional error of 1CP 6 in 
all levels except the highest. The user specifies a 
signal-to-noise ratio S; the MC noise is reduced 
by increasing the number of photons such that the 
fractional error of the levels between iterations is 
smaller than 1/5. For the test problem, 5=100 
was used. 

Complete Linearization (MULTI-type): 
— J. Yates (Rawlings & Yates 2001); In the problems 
presented, the SMMOL code used 400 rays through 
the cloud to compute the intensities at each ra- 
dial grid-point in the spherical cloud using a finite 
difference method to compute the intensities. The 
grid sampling along each ray can adapt simply to 
take account of large changes that affect the op- 
tical depth, typically caused by velocity changes 
along the ray. In this calculation, typically 100 grid- 
points were used to sample the full absorption pro- 
file. The SMMOL codes uses a convergence crite- 
rion of 1 x 10~ 4 for both the level populations and 
the mean intensities. 



4. Description of the test models 

4.1. Models la/lb: a simple 2-level molecule 

The first two test problems (problems la/lb) have a 
simple and cleanly defined setup without velocity gradi- 
ents and using a constant temperature and line width. 
Complicating factors such as a different treatment of the 
spatial grid and frequency sampling are thus kept to a 
minimum, so that every code should in principle be able 
to solve these problems down to their specified accuracy. 

The simple 2-level molecule test setup consists of a 
spherically symmetric cloud with a power law hydrogen 



number density specified by 
nn 2 (r) = n H2 (r ) 



ro 



(15) 



where r is the radius in cm. We take n>H 2 (To) = 2.0 x 
10 7 cm- 3 , r = 1.0 x 10 15 cm and a = -2. The ki- 
netic temperature of this problem is taken to be constant: 
TkinM = 20 K. The abundance of the molecule (which 
we will specify below) is constant as well: X mo \{r) = 
n mo i(r)/nn 2 (r) — X mo i, and the systematic velocity is 
zero everywhere. The spherically symmetric cloud ranges 
from r in = 1.0 x 10 15 cm to r out = 7.8 x 10 18 cm. For 
r < r in and for r > r out the density is assumed to be 
zero. The incoming radiation at the outer boundary is the 
T = 2.728 K microwave background radiation. 

We choose a Active 2-level molecule which is specified 

by 

E 2 -E x =6.0 cm" 1 = 8.63244 K (16) 
92/91 = 3.0 (17) 
A21 = 1.0 x 10~ 4 s" 1 (18) 
K 2 i = 2.0 x 10~ 10 cm 3 s" 1 (19) 

in which the downward collision rate is C21 = ^h 2 ^21 s_1 - 
The total (thermal+turbulent) line width a is given by 
a = 0.150 kms -1 . The line profile is assumed to be a 
Gaussian: 



<f>(y) 



I c 2 (v-v ) 2 
■. exp 1 



ai/oV 71 " 



,2„2 



(20) 



where c is the speed of light in units of km s _1 , and v is 
the frequency at line center. 

We solve the problem for two different abundances: 

Problem la: X mol = 1.0 x 10~ 8 
Problem lb: X mol = 1.0 x 10~ 6 

Problem la is a simple case with a moderate optical depth 
(r ~ 60 from star to infinity at line center). Problem lb 
has much higher optical depth (r ~ 4800) and is therefore 
numerically stiffer. 

This problem was originally described and worked out 
by Dullemond & Ossenkopf (see Dullemond 1999). 
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4.2. Models 2a/2b: A collapsing cloud in HCO + 

The second two problems (problems 2a/2b) are exam- 
ples of problems typically encountered in the field of 
sub-millimeter molecular line modeling. These problems 
have much of the complicated physics included, such as 
velocity- and temperature gradients, non-constant line 
widths, multiple levels etc. Moreover, the problem is based 
on a numerical model provided externally and there- 
fore introduces the complication of a given intermediate- 
resolution spatial grid. The resulting spread of the results 
therefore gives a good indication of how accurate the pre- 
dicted line strengths and shapes for a typical every-day-life 
model calculation are. 

The model on which these "realistic" test problems 
are based was constructed by Rawlings et al. (1992, 1999) 
to analyze HCO + data for an infalling envelope around 
a protostar. The prototypical example of such a 'class 0' 
young stellar object is B335. The model describes how the 
cloud core collapses from the inside-out. Starting from a 
nearly isothermal sphere in pressure balance, a perturba- 
tion triggers the center of the cloud to collapse. This sends 
out a rarefaction wave at the local speed of sound, which 
causes the outer parts of the cloud to collapse as well. 
The model is similar to the analytical inside-out collapse 
model by 3hu (1977), but includes more realistic physics. 
No rotation is assumed, which makes it ideal for a 1-D 
spherical comparison. 

Figure || shows the structure of the cloud at a par- 
ticular time during the collapse phase. This is the input 
model for our test cases. The collapse can clearly be seen 
in the radial velocity which is for radii greater than 
10 17 cm, and directed towards the source for smaller radii. 
The density profile is given by a power-law of the form 
n(r) = no(r/ro) m , where m = —1.5 inside the collapsing 
sphere (r < 10 17 cm) and m = —2.0 outside. The model 
parameters are specified at 50 radii logarithmically spaced 
between 10 16 and 4.6 x 10 17 cm. 

The infall of a cloud has a significant influence on the 
molecular line emission. The lines will be skewed to the 
blue, or double peaked with a stronger blue-shifted com- 
ponent, depending on the velocity field and optical depth 
of the line (e.g. Myers et al. 2000). This feature can be 
intuitively understood because the foreground gas has a 
lower excitation temperature than the background gas. 
This only holds when the foreground gas is optically thick 
enough in the line. Otherwise a single-peaked lineshape, 
centered on the source velocity, results. 

For the comparison presented here, the ion HCO + is 
used, where an abundance compared to H2 is adopted 
of lxl0~ 9 for problem 2a, and lxl0~ 8 for problem 
2b. This causes the optical depths r for the lowest 4 
transitions to range from 0.1 to 10 for problem 2a, 
and from 1 to 100 for problem 2b. A common file of 
Einstein- A coefficients and collisional rate coefficients of 
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Fig. 2. Physical structure of the test model 2. Upper 
left: density lower left: velocity upper right: temperature, 
lower right: turbulent line width b [km s -1 ]. Changes in 
the density and velocity distributions are seen at the point 
where infall starts. The temperature of the gas at the in- 
side rises due to the infall while the outside is heated by 
the interstellar radiation. 

of the HCO + ion are included in the model. Since the 
rate coefficients are temperature dependent, all downward 
(upper— dower) rate coefficients are given for a number of 
temperatures (10, 20, 30, 40 K). At each temperature, 
a linear interpolation is performed to calculate the lo- 
cal downward collisional rate coefficients. Performing this 
interpolation accurately is extremely important because 
the exponential relation between the excitation and de- 
excitation coefficients can introduce large deviations. The 
upward (lower— >upper) collisional rate coefficient is sub- 
sequently calculated using the detailed balance relation, 



_ 9i c E ul /kT 



(21) 



HCO+ with H 2 in J=0 flMontciro 1985| , |Green 1975[) has 
been used in all models. These can be downloa ded from 
http : / / www . strw . leidenuniv . nl/ ^radtrans] . 21 levels 



The fact that the molecular data, in particular the 
collisional rate coefficients, are not known with infinite 
precision introduces an additional uncertainty in the so- 
lution of the radiative transfer equations. The use of 
different collisional rate coefficients has an impact on 
the level populations and thereby on the strength of 
the predicted emission lines. However, the spread in 
available rate coefficients in the literature has no ef- 
fect on the comparison presented here, as all the codes 
use the same molecular data as input. On a webpage 
tittp : //www, strw. leidenuniv .nl/^moldata molecular 
files for a large number of molecular species will be made 
available in the future (Schoier et al., in preparation), in- 
cluding a comparison of different rates from the literature 
and their effect on the model results. 
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With these sets of problems the different codes are 
tested on the following features: 

— Convergence of the code 

— Radiative transfer in optically thin lines 

— Radiative transfer in optically thick lines 

— Sampling of the radiation transfer in the presence of a 
velocity field 

— Incorporation of the cosmic microwave background ra- 
diation (Tcmb=2.728 K), which influences the lower 
levels near the outside of the cloud. 



5. Results 

5.1. Results of problems la/lb 

Test problems la and lb are relatively simple in the sense 
that no velocity gradients and temperature variations are 
present. Also, since the problem setup is specified analyt- 
ically, the grid resolution can be chosen to be as high as 
required for the particular code the participant is using. 
By virtue of the zero systematic velocity, the treatment 
of the inner boundary is not much of an issue, because 
the line is perfectly thermalized at that radius. These test 
problems therefore purely test the radiative transfer as 
such. 

The results for the problems la and lb are shown in 
Fig. ||. The interior of the cloud is completely thermalized 
(T ex = Tkin)- The decline of 77,2 at larger radii is due to 
the departure from LTE. At the outer radii the population 
again saturates to a constant value, which is due to the 
microwave background radiation (T cx = Tomb)- All codes 
agree reasonably well, although small differences remain 
distinguishable. The relative differences between the codes 
are shown in Fig.|| in comparison to the mean. 

If we consider the 'mean' solution to be the correct 
one, then the typical errors are of the order of a few per- 
cent or less. The errors are clearly greater for problem lb 
than for problem la. This is to be expected, since prob- 
lem lb has higher optical depth and is therefore numer- 
ically more challenging. For problem la, the codes show 
a small spread with a standard deviation of 0.2 %. i.e. 
they agree with each other to within 2 • 10~ 3 relative dif- 
ference, comparable to the accuracy criteria of the codes, 
clearly visible in the noise for the Monte Carlo results. For 
problem lb the best agreement reached between the codes 
shows a standard deviation with a local maximum of 2% 
(An 2 /n 2 ~ 2 ■ 1(T 2 ). To achieve maximum accuracy, the 
codes used a spatial resolution of (Ri+i — Ri)/Ri ~ 0.05 
for this problem, but this high spatial resolution did not 
reduce the spread between the codes to the 10 -3 level. 

For problem lb (the most difficult of the two) an 
ALI code with Ng acceleration with spatial resolution 
(Ri+i — Ri)/Ri ~ 0.04 typically takes about 160 itera- 
tions to meet the converge criterion of An2/n,2 < 10~ 7 . 
The same problem, but for a less stringent convergence cri- 
terion (Att^/tt^ < 10~ 4 ) and for a lower spatial resolution 
— Ri) I Ri — 0.2) typically takes about 40 iterations. 




1 b 




Fig. 3. Results of the two-level molecule problems. Shown 
here are the populations of the upper level of the 2-level 
molecule. Top panel: problem la. Bottom panel: problem 
lb. Plotted are the (A)MC codes, denoted by solid lines 
and the ALI/SMULTI codes using dashed lines. 



The relative error introduced by the lower spatial resolu- 
tion and lower convergence criterion for this problem is 
at most 1%, provided that second-order (or higher) inte- 
gration of the transfer equation is used. The use of first- 
order integration introduces errors of the order of 12%, 
even for the high-resolution runs. For test problem la, 
at low resolution and with a weak convergence criterion 
(Ari2/ri2 < 10~ 4 ), an ALI code typically uses 18 itera- 
tions. 

For the AMC-type codes with Ng acceleration (Juvela 
1997), problem lb typically requires about 75 iterations, 
with about 120 photon packages (rays) per iteration. It 
was found that the application of Ng acceleration at the 
very start of the iteration had a negative influence on the 
convergence rate. Good convergence was achieved by ap- 
plying Ng only after about 25 iterations. As in the case of 
the ALI type codes, problem la was much more benign: 
15 iterations with about 120 photon packages (rays) per 
iteration. 

It should be said that the AMC-type codes allow con- 
siderable freedom in how to disentangle the number of 
iterations and the number of photon packages per cell per 
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Fig. 4. As Fig.g, but now shown as the relative difference from the average of the results for problems la (left) and 
lb (right). Plotted are the (A)MC codes, denoted by solid lines and the ALI/SMULTI codes using dashed lines. The 
notation for each of the codes is given in Table 1. Note that the differences between the codes is of the same order as 
the Monte Carlo noise for problem la, indicating that the codes have converged to the correct level. 



iteration. Therefore the numbers given here may vary from 
code to code, i.e., in most codes the number of photon 
packages per iteration per cell is used. For the accelerated 
Monte Carlo codes, acceleration of the convergence only 
operates when some of the cells are optically thick. If all 
cells are optically thin but the entire source optically thick, 
the local radiation field becomes negligible compared to 
the overall radiation field and separating local and global 
contributions no longer effectively speeds up convergences. 

The SMULTI code converged to below 1CT 4 in 30 it- 
erations and to below 10~ 5 in 50 iterations for the second 
model. For the first model the code converged in 10 iter- 
ations. 

5.2. Results of problems 2a /2b 
5.2.1. Level populations 

The results for the test problems 2a and 2b are shown in 
Fig. |. The direct comparison is shown for two different 
levels of HCO + , J=l and J=4. These levels were chosen 
because they represent the emitting levels of two easily 
and often observed lines of the molecule probing different 
density regimes. The J=l-0 line lies at millimeter wave- 
lengths at 89 GHz and the critical density of the J=l 
level is 3.4 x 10 4 cm -3 . The J=4-3 line occurs in the sub- 
millimcter at 356 GHz, with a critical density for exciting 
the J=4 level of 1.8 x 10 6 cm~ 3 . The different model re- 
sults are plotted on top of each other, together with one 
additional calculation where the CMB radiation field was 
deliberately ignored. In the level populations, the effect 
of the CMB is visible as a lack of excitation of the J=l 
level in the outer regions. The figure immediately shows 
that in the outer region, the excitation at line center of 
the J=l-0 line is controlled by the CMB field and not 
by the local temperature and density. The J=l level lies 
only 4.3 K above ground, close to the Tomb (2.728 K) 
temperature, so that the 1-0 transition can be effectively 



excited by the peak of the CMB radiation. The density in 
the outer regions of the collapsing cloud is too low for the 
molecule to be excited to LTE. The J— A level population 
is less affected by the CMB radiation field. 

To show a more quantitative measure of the accuracy 
of the results, the level populations are plotted versus their 
difference to the mean of all the results (Fig. ^). Only re- 
sults including the CMB radiation field are used to cal- 
culate the mean. At the inner boundary, the solutions di- 
verge into two main groups. This is a result from the in- 
ner boundary condition adopted by the different authors. 
The density and temperature of the test problem were not 
specified within the inner radius and were either chosen 
as an empty sphere or a non-rotating sphere of constant 
density. The different solutions near the inner boundary 
have a relatively small error (Fig. ||) showing that in this 
case the inner boundary has little influence on the over- 
all solution and can be ignored in the specification of the 
problem. 

The calculated populations of the J=l level show a 
standard deviation ranging from 2% close to the star down 
to 1.5% for problem 2b and below 1% in problem 2a. The 
J=4 level has a larger spread with a standard deviation 
of 8% close to the inner boundary up to 12% at r=2x 10 17 
cm and down to 2% at the outer boundary for the high 
t (problem 2b) case. The HCO + lines have high critical 
densities and the J=4 level population is far from LTE. 
In addition, the optical depth is high, making the calcula- 
tion for the 4-3 line particularly difficult. In problem 2a 
the J = 4 level rises from 3% close to the inner bound- 
ary up to 4% in the transition region and down to 1% 
close to the outer boundary. Lacking an analytical solu- 
tion, we take the average of the numerical results as the 
best estimate of the exact solution. The deviations from 
the average therefore give an estimate of the implicit and 
explicit approximations of the codes. Fitting models to ob- 
servational data with criteria better than these deviations 
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Fig. 6. Differences of the populations in the J=l and J=4 levels to the mean for both model 2a and 2b. The solid 
lines represents the (A)MC based codes and the dashed lines the ALI based codes. The notation of each of the lines 
is given in Table 1. Only codes A-G participated in this comparison. 



will not lead to more accurate estimates of the physical 
parameters. 

For all the problems, the level populations were also 
calculated using LI by Dullemond and Wiesemeyer. The 
populations agreed within the errors as shown in Figs. 4 
& 6 for the low r cases. However in both the lb and 2b 
problems the solution was not reached using similar con- 
vergence criteria as the ALI calculations. Even though the 
LI method itself should in principle find the same solution, 
the number of iterations and more strict convergence cri- 
teria make it a numerically too costly task to solve. 

The problem is particularly sensitive to the gridding of 
the physical parameters. A comparison of two runs, where 
one had twice as many gridcells, showed a significant de- 
crease of the error compared to the mean value for one of 
the Monte Carlo code (C). This shows the importance of 
a correct gridding of the problem. 

5.2.2. Excitation temperature 

As a second comparison, the excitation temperature 
(Fig. [7]) is plotted for all solutions, which may be a more 
intuitive manner to show the results. When level popu- 
lations are in LTE, the excitation temperature equals the 
local kinetic temperature. Comparing the results with Fig. 



, the excitation temperature of the 1-0 line is close to 
LTE near the inner boundary but becomes subthermal 
at larger radii. The J=4-3 excitation temperature is al- 
ways well below its LTE value, which is due to its much 
higher critical density. As Fig. [?] shows, the excitation tem- 
perature never drops below 2.728 K, except in the case 
where the CMB radiation was ignored. In the inner re- 
gion (r < 5 x 10 16 cm), the excitation of the J=l level is 
dominated by collisions, while in the outer parts, it is dom- 
inated by radiation. The transition region between these 
two extremes shows the largest differences in the results 
(Fig. la). 

5.2.3. Line profiles 

Figure || compares the results in terms of the line pro- 
files. The calculated level populations at each position in 
the cloud are used to compute the velocity profiles of the 
selected lines using a program which calculates the sky 
brightness distribution. To ensure an equivalent emitting 
mass, the inner sphere is assumed to be empty for all mod- 
els. The profiles are calculated by constructing a plane 
through the origin of the cloud perpendicular to the line 
of sight, with a spatial resolution small enough to sample 
the physical distributions. A ray-tracing calculation is per- 
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Fig. 5. Populations of levels J=l and J=4 in the optically 
thin (model 2a) and optically thick (model 2b) cases. The 
dotted line is the solution where no cosmic background 
radiation was added, the solid lines represent the (A)MC- 
based codes and the dashed lines the ALI-based codes. 
Due to the low temperatures and densities of the problem, 
most molecules are in the lowest rotational states. Only 
codes A-G participated in this comparison. 




W 
10 

g 25 

8 20 



8.15 
g 



10 



3 5 

u 
x 

10 



16 


10 17 10 




2b 




J=1 /J=0 







16 



17 



10 

r fcrnl 



10 




r Icml 



Fig. 7. Excitation temperatures [K], as defined in Eq. (g), 
for the levels J=1/J=0 and J=4/J=3. Top panels: low 
optical depth case (problem 2a); bottom panels: high opti- 
cal depth case (problem 2b) . The dotted line indicates the 
model result with the CMB radiation neglected, the solid 
lines represent the (A)MC-bascd codes and the dashed 
lines the ALI-based codes. When the lines arc in LTE, the 
excitation temperature is equal to the kinetic temperature 
(Figure |). 



formed through this plane from — oo to +oo, keeping track 
of the intensity in the different velocity bins. For this cal- 
culat io n the program SKY was used, part of the RATRAN 
code ( jittp : //talisker . as . ar izona . edu/ ^michiel/ ) 
By applying a single common code for these calculations, 
no secondary uncertainties are introduced. The resulting 
sky brightness distribution is convolved with a beam of 
14" for the J=4-3 transition, appropriate for the James 
Clerk Maxwell Telescope (JCMT) at this frequency. For 
the J =1-0 transition, the beam size is taken to be 29" 
representative of the IRAM 30m telescope. The distance 
is assumed to be 250 pc. 

Comparison of Fig. ||a and c shows that increasing the 
abundance of HCO + by a factor of 10 changes the 3=1- 
peak intensity by only a factor of 2.5. The J =1-0 line 
becomes more self-absorbed, indicating that the line has 
indeed become more optically thick. The J=4-3 emission 
line is optically thin for the low abundance, as seen from its 
single-peaked profile in Fig. ^Jd. The more optically thick 
model shows the characteristic asymmetric line structure, 
skewed to the blue. The differences in the level populations 
are visible in the line profiles. In the J =1-0 line profiles 
(Fig. ||c) all solutions lie on top of each other in the line 
wings and only a slight difference can be seen at line cen- 
ter. The error is largest in the region with low velocities 
represented by the larger errors in the center of the line 
profile. The integrated intensity profiles differ by a few % 



for low t and 6-7% for the high r case. The error is in gen- 
eral larger for the J=4 level, which shows up as a larger 
error in the integrated lines for the J=4-3 transition. For 
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Fig. 8. Calculated intensity [K] of the J =1-0 and J=4-3 
lines for the problems 2a and 2b. In the high optical depth 
case (lower two panels), the blue-shift or 'infall asymme- 
try' is visible. The low resolution of the spatial grid be- 
comes apparent in the line emission at ± 0.75 km s _1 . 
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low r the error is w 2%, but for high optical depths de- 
viations up to 12% are found for the most outlying case. 
Although substantial, these differences are generally less 
than the current calibration uncertainties of the observa- 
tional data, which are typically 20% to 30%. 

The differences in the profiles shown here give a rough 
indication of the quality of fitting necessary to interpret 
observational data. It should be noted that some differ- 
ences are likely due to different implementations of the 
cloud model. One can fit a model precisely to an observa- 
tional data set using a particular radiative transfer code; 
however the results will not necessarily be the same if a 
different code is used. The derivation of physical param- 
eters is therefore always limited in accuracy to the error 
bars given by the code itself in comparison to other codes. 

6. Discussion and conclusion 

In spite of the fact that the spread in our solutions lies 
within the accuracy limits of current-day (sub-)millimeter 
instruments, it is important to understand how this spread 
comes about. The first two test cases were defined analyt- 
ically, and avoided difficulties of gridding as much as pos- 
sible. The small spread of these results, even at high op- 
tical depths, seems to show that in principle the radiative 
transfer is done correctly by all codes. The larger spread 
in the HCO + collapsing cloud problem can therefore be 
attributed to the problem of gridding, both in space and 
in frequency. The coarse spatial gridding in the setup of 
the problem, and the presence of strong velocity gradi- 
ents makes it very likely that different interpretations of 
the sub-grid behavior of temperature, density and velocity 
lead to slightly different results. 

The main results of the comparison of radiative trans- 
fer codes for molecular lines can be summarized as follows. 

— All relevant methods currently used in molecular as- 
trophysics agree to a few % for optical depths up to 
r ^4800. At high optical depth and in transition layers 
between collision-dominated and radiation-dominated 
excitation, relative differences up to 2 % can arise for 
well defined "simple" problems. In practice, for more 
complex models, the relative differences are higher due 
to gridding and geometry problems, reaching relative 
differences up to 12 % locally for the high optical 
depth model described in this paper (problem 2b) . 
However, in most practical applications in molecular 
astrophysics the optical depths are lower than in our 
test problems, and therefore this 12% error can be re- 
garded as an upper limit. 

— The error bars on the line profiles are generally much 
less than 10%, and therefore lie within the calibration 
errors of typical (sub-)mm observations. Only in the 
worst case these errors may be comparable to the ob- 
servational uncertainty. 

— The choice of gridding is of extreme importance, and is 
one of the major causes of the deviations in the prob- 



lems presented. Special care should be taken when con- 
structing models for specific problems. 

— Abundances and excitation temperatures derived from 
lines which are formed in the transition layer should 
be interpreted with caution. 

— The direct comparison of results of different programs 
speeds up significantly the debugging process of new 
programs. 

The work presented here is one step in the direction 
of standardization of radiative transfer computations in 
molecular rotational lines. Suggestions for future work 
in this area include: (1) establishment of a database for 
collisional rate coefficients, accessible to the community 
through the WWW; (2) a campaign to have inelastic rate 
coefficients measured or calculated for those molecules for 
which these data are unavailable, and (3) comparison of 
2D radiative transfer codes, with a test problem based on, 
e.g., the rotating flattened collapse described by Terebey 
et al. (1984) or the "sheet" models of protostellar collapse 
by Hartmann et al. (1996). These developments will be 
vital to interpreting the high-quality data which e.g. HIFI 
(the heterodyne instrument onboard the Herschel Space 
Observatory) and ALMA will provide. 
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